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Abstract 

We investigate a projective integration scheme for a kinetic equation in the limit 
of vanishing mean free path, in which the kinetic description approaches a diffusion 
phenomenon. The scheme first takes a few small steps with a simple, explicit method, 
such as a spatial centered flux/forward Euler time integration, and subsequently projects 
the results forward in time over a large time step on the diffusion time scale. We show 
that, with an appropriate choice of the inner step size, the time-step restriction on the 
outer time step is similar to the stability condition for the diffusion equation, whereas 
the required number of inner steps does not depend on the mean free path. We also 
provide a consistency result. The presented method is asymptotic-preserving, in the 
sense that the method converges to a standard finite volume scheme for the diffusion 
equation in the limit of vanishing mean free path. The analysis is illustrated with 
numerical results, and wc present an application to the Su-Olson test. 



1 Introduction 



In many applications, ranging from radiative transfer over rarefied gas dynamics to cell 
motion in biology, the underlying physical system consists of a large number of moving and 
colliding particles. Such systems can be accurately modelled using a kinetic mesoscopic 
description that governs the evolution of the particle distribution in position-velocity phase 
space. In a diffusive scaling, when the mean free path of the particles is small with respect 
to the (macroscopic) length scale of interest, a macroscopic description involving only a 
few low-order moments of the particle distribution (such as a diffusion equation in neutron 
transport and radiative transfer, or fluid equations in rarefied gas dynamics) may give a 
rough idea of the behavior. However, refining the description is uneasy because, whereas the 
physical model becomes much simpler in the diffusion limit, a direct numerical simulation 
of the kinetic model tends to be prohibitively expensive due to the additional dimensions 
in velocity space and stability restrictions that depend singularly on the mean free path. 
We consider dimensionless kinetic equations of the type 

dtf, + --V,f, = ^Q{f,), (1.1) 

modeling the evolution of a particle distribution function f^{x,v,t) that gives the distri- 
bution density of particles at a given position x U C M.'^ with velocity v ^ V C M.'^, 
d > 1, at time t, the collisions being embodied in the operator Q. The parameter e > 
is meant as the ratio of the mean free path over the characteristic length of observation, 
i.e. the average distance traveled by the particles between collisions. The diffusion limit is 
obtained by taking e — t- 0. Under some appropriate assumptions, which will be detailed 
in section 2, the unknown then relaxes on short time-scales to an equilibrium, in which 
the dependence on v is fixed, and the dynamics of the system on long time-scales can be 
described as a function of the density pe{x,t) = {fe{x,v,t)), where 

(•>= / -dpiv), 
Jv 

is the averaging operator over velocity space and {V, dp) denotes the measured velocity 
space. For e — )• 0, the density p = lim^-i-o Pt satisfies formally the diffusion equation 

dtp-d/\^p = Q, dp = {v^). (1.2) 

The difficulty in studying the asymptotic behavior numerically is precisely due to the 
presence of two time-scales. On the one hand, explicit time integration of equation (1.1) is 
numerically challenging because one is forced to take very small time-steps 5t when e — )■ 
to stably integrate the fast relaxation. Indeed, due to stability considerations, 6t needs to 
be shrunk as e — )■ to properly satisfy both the e-dependent hyperbolic Courant-Friedrichs- 
Lewy (CFL) condition for equation (1.1) and the stability constraints for the collision term. 
On the other hand, implicit schemes are computationally expensive because of the extra 
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dimensions in velocity space. At the same time, the equation closely resembles a diffusion 
equation in that limit, for which a parabolic CFL condition of the type At = O(Ax^) 
(independently of e) would be desirable. 

A number of specialized methods that are asymptotic-preserving in the sense introduced 
by Jin [11] have been developed that can integrate equation (1.1) in the limit of e — t- with 
time-steps that are only limited by the stability constraints of the diffusion limit. We briefly 
review here some efforts, and refer to the cited references for more details. In [13, 16], 
separating the distribution / into its odd and even parts in the velocity variable results 
in a coupled system of transport equations where the stiffness appears only in the source 
term, allowing to use a time-splitting technique [24] with implicit treatment of the source 
term; see also related work in [11, 12, 17, 18, 22]. When the collision operator allows for an 
explicit computation, an explicit scheme can be obtained by splitting / into its mean value 
and the first-order fluctuations in a Hilbert expansion form [9] under a classical diffusion 
CFL. Also, closure by moments [5, 21, e.g.] can lead to reduced systems for which time- 
splitting provides new classes of schemes [4]. Alternatively, a micro- macro decomposition 
based on a Chapman-Enskog expansion has been proposed [20], leading to a system of 
transport equations that allows to design a semi-implicit scheme without time splitting. 
An innovative non-local procedure based on the quadrature of kernels obtained through 
pseudo-differential calculus was proposed in [2]. 

Our goal is to introduce a different point of view, based on methods that were devel- 
oped for large multiscale systems of ODEs. In [8], projective integration methods were 
introduced as a class of explicit methods for the solution of stiff systems of ordinary differ- 
ential equations (ODEs) that have a large gap between their fast and slow time scales; these 
methods fit within recent efforts to systematically explore numerical methods for multiscale 
simulation [6, 7, 14, 15]. In projective integration, the fast modes, that correspond to the 
Jacobian eigenvalues with large negative real parts, decay quickly, whereas the slow modes 
correspond to eigenvalues of smaller magnitude and are the solution components of practical 
interest. Such problems are called stiff, and a standard explicit method requires time steps 
that are very small compared to the slow time scales, just to ensure a stable integration of 
the fast modes. Projective integration circumvents this problem. The method first takes a 
few small (inner) steps with a simple, explicit method, until the transients corresponding to 
the fast modes have died out, and subsequently projects (extrapolates) the solution forward 
in time over a large (outer) time step; a schematic representation of the scheme is given 
in figure 1. A stability analysis in the ODE setting was presented for a first order version, 
called projective forward Euler [8], and an accuracy analysis was given in [26]. Higher-order 
versions have been proposed in [19, 23]. 

Projective integration methods can offer a number of important advantages for the 
simulation of kinetic equations. In particular, they are fully explicit and do not require 
any splitting, neither in time, nor in microscopic and macroscopic variables. In this work, 
we will analyze the properties of projective integration for kinetic equations on diffusive 
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Figure 1: Schematic representation of projective forward Euler. In a stiff system with a spectral 
gap, the values of the solution obtained by straightforward explicit integration (dots) are quickly 
attracted to a slow manifold (solid curve). Projective integration takes a few explicit inner steps 
until the solution has come down to this slow manifold. Then, a chord slope estimation is performed, 
and a big projective outer step is taken. Since the result of the outer step does not lie on the slow 
manifold, a number of inner steps is taken again, and the procedure is repeated. 

time scales in one space dimension, keeping in mind that these methods extend readily to 
higher dimensions. The paper is organized as follows. In section 2, we give the necessary 
preliminaries on the model problem (1.1) and present its diffusion asymptotics. We also 
introduce the Su-Olson test case that will be used in the numerical experiments. In section 
3, we develop the numerical scheme. We present the brute-force inner schemes and the 
projective outer forward Euler method. In section 4, we show that, when choosing the 
inner time step 6t = e^, the stability condition on the outer time step is independent of 
e, and similar to the CFL condition of the limiting heat equation. Moreover, the required 
number of inner steps is also independent of e when e — )■ 0. Subsequently, in section 5, we 
study the consistency of the method and, using the stability condition derived in section 4, 
we give a bound of the convergence error, enabling to conclude that the presented scheme 
is asymptotic-preserving. We provide numerical illustrations for the linear model and the 
Su-Olson test in section 6. Finally, section 7 contains conclusions and an outlook to future 
work. 

2 Model problem and diffusion asymptotics 
2.1 Linear relcixation 

We consider equation (1.1) in one space dimension, 



that can be interpreted as the difference of a gain and a loss term. In the remainder of 
the text, we restrict the discussion to a periodic setting in one space dimension, hence 



dtfe + -d.fe = ^Qife) 




and specify the collision operator as a linear relaxation 



= (/e> - /e = /9e - /. 
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X £ T = [0,1) and v G y C M. The results of the derivation, however, could easily be 
generalized to higher space dimensions. 

Throughout the paper, the measured velocity space (V, fi) is required to satisfy 

jydfiiv) = 1, 

jy h{v)dfi{v) = 0, for any odd integrable function h : V — > M, 
Jy v'^dfi{v) = d> 0. 

Typical examples are 

• V = (—1,1) endowed with the normalized Lebesgue measure, for which we have 
d = 1/3; 

• the discrete velocity space 

2j — 1 

V = {-Vp, . . . , -vi,vi, . . .,Vp}, with Vj = — , j e J = J_^_u J-, (2.2) 

/p 

where J'± := {±1, . . . , =bp}, endowed with the normalized discrete velocity measure 
1 

dfJ'iv) = — J2jeJ "^(^ ~ ^j)' which 
2p 



^ 2p 12p2 



so that dp — )• 1 /3 as p — )■ oo; 



• y = M endowed with the Gaussian measure d^{v) = {2-k) exp(— ?;^/2)ci?;, for which 
d=l. 

These assumptions on the velocity space result in a number of properties for Q, namely : 

1. Q is a bounded operator on L'p{V, dfi{v)), I < p < oo; 

2. Q is conservative, i.e. 

\ffeLHv,df,iv)):{Qif))=0; 

3. the elements of the kernel of Q are independent of v. 

As e — 7- 0, the frequency of collisions increases; hence, we may formally propose to write 
fe as a perturbation of the macroscopic density = (f^) using a Hilbert expansion: 

fe{x,v,t) = pe{x,t) + ege{x,v,t). (2.3) 

Equation (2.1) can then alternatively be written as 

dtPeix,t) + {vdxgeix,v,t)) = 0, 

1 1 (2.4) 

dtge{x,v,t) + - {vdxge - {vd^ge)) = — o (^e + vdxPe) ■ 
e 
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Taking formally the limit e — )■ yields 

dtp,-{v^)d^^p, = 0{^), (2.5) 

since 

= -vdxPe + e{{vd^ge) - vd^ge} + 0{e^). (2.6) 

The approximation (2.5) is consistent at order 2 in e with the diffusion equation (1.2). We 
refer to [5] for a recent starting point of the literature on the convergence of /e to p, solution 
of (1.2). 

For concreteness, we will use the discrete velocity space (2.2) in our analysis and sim- 
ulations. Thus, we are now interested in a vector- valued function : x T — > M?^, of 
which we denote the component corresponding to Vj as fej{t,x), Vj S J^. In this setting, 
the density is given as pe = (YIj f<^,j^ /C^P)- The kinetic equation (2.1) then reads 

VjGj, dtU, + ^d.U, = ^^^/^. (2.7) 
2.2 Su-Olson equation 

While the linear kinetic equation (2.1) is an ideal model problem for analysis purposes, we 
will also show numerical results for a more challenging test case, namely the traditional Su- 
Olson benchmark, a prototype model for radiative transfer problems. Here, the unknown 
/e represents the specific intensity of radiations, which interact with matter through energy 
exchanges; see e.g. [2, 3, 4, 25]. A complete model couples a kinetic equation for the evolution 
of /e with the Euler system describing the evolution of the matter. In the Su-Olson test, this 
coupling is replaced by a simple ODE describing the evolution of the material temperature. 
The system reads 

dtfe + -fe = \ {Pe " fe) + (6 - p,) + S, 

e (2.8) 

dtO = aa {pe - e) . 

Here, O = T^, with T the material temperature, and S" is a given source depending on x. 
In our simulations, the parameter (Tq = 1. 



3 Numerical scheme 

3.1 Finite volume formulation 

We consider a uniform, constant in time, periodic spatial mesh with spacing Ax, consisting 
of cells Ci = 1 < i < Nx with NxAx = 1 centered in Xj, where Xi = iAx, 

and a uniform mesh in time = [t^,t''~^^), k > = k6t where 6t is the time step. We 
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adopt a finite volume approach, that is we integrate (2.7) on a cell Mj^^ = CixTj. to obtain, 
Vj G J, 

(^,■(^^•+^•)-/e,,(t^•)) + ^f / (/.,,(-,Xm/2)-/ej(-,X,_i/2))) =4 / (Pe-/e,,). 



c, 

In order to simplify the notations, the solution of the continuous equation (2.7) will 
always be denoted with the subscript e whereas the solution of the numerical scheme will 
be denoted / = [fij]- This leads to the conservative forward Euler scheme 

=fh-7L (^(^)ti/2. - ^(/)ti/2.) + § (pf - ft^ , (3.1) 

where ff^ denotes an approximation of the mean value fej{t'' , x)dx, the numerical flux 
<i^{f)'lj^il2 j is an approximation of the flux at the interface 2:4+1/2 a-t time i*^ in the equation 
for the velocity j, and = (/f ), the average being taken over the velocity index. We will 
consider upwind or centered numerical fluxes that are given by 



^ I Vjff,, if ?■ G J+ (vj > 0), 

Mf)li/2, = ' / (3.2) 



respectively. We introduce a generic short-hand notation, 
with / G M^^^^p ^ square matrix of order Nx x 2p. 

Remark 3.1 (Maximum principle). We recall that, while the centered flux scheme (3.3) 
is second-order accurate in space, it does not obey a maximum principle, and hence may 
lead to unphysical oscillations, the centered transport scheme being violently unstable for 
transport equations. Any projective integration scheme based on the centered flux scheme 
should therefore also violate the maximum principle. However, since the kinetic equation 
(2.1) is consistent at order 2 in e with a diffusion equation (2.5), the oscillations are quickly 
stabilized as e — )• (see also the discussion on consistency in section 5). 



3.2 Projective integration 

Because of the presence of the small parameter e, the time steps that one can take with the 
upwind scheme are at most 0(e), due to the CFL stability condition for the transport equa- 
tion, or O(e^) due to the relaxation term. However, in the diffusion limit, as e goes to 0, the 
equation tends to the diffusion equation (1.2), for which a standard finite volume/forward 
Euler method only needs to satisfy a stability restriction of the form At < Ax^/(2(i). 
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In this paper, we consider the use of projective integration [8] to accelerate brute force 
integration; the idea is the following (see also figure 1). Starting from an approximate 
solution at time = NAt, one first takes K + 1 inner steps of size St, 

fNMi^Sstf'''', k = 0,...,K, 

in which the superscript pair (A^, k) represents an approximation to the solution at t^''^ = 
NAt + k5t. The aim is to obtain a discrete derivative to be used in the outer step to 
compute f^~^^ = f^+^'^ via extrapolation in time, e.g., 

fN,K+l _ fN,K 

jN+l ^ fN,K+l ^ -(K+ l)6t) ^ . (3.5) 

ot 

This method is called projective forward Euler, and it is the simplest instantiation of this 
class of integration methods [8]. Adams-Bashforth or Runge-Kutta extensions of (3.5), 
giving a higher order consistency in terms of At, are possible [19, 23]. 

Projective integration is a viable asymptotic-preserving scheme if, as e goes to 0, we 
have (i) a stability for the outer time step At that should satisfy a condition similar to the 
CFL condition for the diffusion equation, (ii) a number of inner steps that is independent of 
e and (iii) the consistency with the diffusion equation (1.2). The analysis of these properties 
will be performed in the next sections. 



4 Stability analysis 
4.1 Notations 

To perform a Von Neumann analysis of the projective forward Euler scheme, we need the 
following notations : 

• e := (1, . . . , 1)'^ G M^^' and V := ee'^ is the orthogonal projection on Span(e); 

1 

• yW £ R^P, VW = {W){l,...,lY = ^{W)e, with {W) = — ^Xj, 

Slllf ^) 

• for all C = CAx, ^ G M, Vc(C) := diag(7;p, . . . ,vi,-vi, . . . , -Vp) and 

^ 2sin(C/2) ^.^g( . . . ^ -v,e-<l\ -v^e-^l^). 

Ax 

We also introduce the symbol V (a, /?) to denote a closed disk with center a and radius j3. 



4.2 Forward Euler schemes 

Let us first locate the spectrum of the matrix Sst defined in (3.4). We denote h = Sstf and 
compute the Fourier series in space of periodized reconstructions of / and h as constant- 
by-cell functions F : x fi if x £ Ci and H : x hi if x £ Ci : Vm G 



Him) = AF{m) = 
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where F{m) := e~'^'^'^'^^ f{x)dx, and correspondingly H{m), are the m-th Fourier coef- 
ficients of F and H and V is the (diagonal) Fourier matrix of the finite volume operator 
chosen for the convection part. We will give hereafter the results for V = Vc and V = Vu- 
Since F G L^^O, 1) ^ {F{m))m S ^{1.) is an isometry, studying the stability of the scheme 
is equivalent to studying the spectrum of A. 

We first prove an auxiliary result. For the sake of simplicity, the matrix V being a 
diagonal complex matrix, we can study the spectrum of a typical matrix A = D + P where 
P = (1, . . . , 1)^(1, •••,!) and Z) is a diagonal matrix with complex entries. 

Proposition 4.1. Let D = diag(L'i, . . . , D2p), with Dj G C. Then the following properties 
of A = D + P hold : 

(PI) if Dj = Dk implies j = k (HI), the eigenspaces of A are all of dimension 1 and the 
spectrum of A does not contain any D^ ; 

(P2) if moreover we assume that Dj = I?2p-j+i Vj G J+ (H2), then an eigenvalue of A is 

— either real 




(a) (P2) (b) (P3) 



Figure 2: Spectrum of A in case (P2) and in case (P3) where 1)^=1) (O, e max^g j+ (ja^ | + |/5j|))-+ 



(P3) if, in addition to the previous hypotheses, D is of order e, e being small, that is 
D = e diag(ap — i/3p, . . . , qi — i/3i, ai + i/3i, . . . , Op + ifip) (H3), then 

Sp(^) C {v (^0,e inax(|aj| + \fij\)^ \ U {A(e)} 

where the only real eigenvalue A(e) is simple and can he expanded as 

A(6) = 2p (l + 6^ - ^((a - {a)f + P^)^ + o{e'), 

(see figure 2.b). 
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Proof. Assume (HI). Let {X,W) be an eigenvalue and an associated eigenvector. Then 
{D + P)W = XW {Xhp - D)W = yfIp{W)e. There are two possibihties : 

• either {W) = and there exists a pair of indices (/ci, ^2) G {1, . • • , 2p}, ki 7^ ^2 such 
that Xfc^ 7^ and Xfcj 7^ 0, which imphes that A = D^-^ = Dk2, which is incompatible 
with (HI) ; 

• or (W) / and W = ^/2p{W){XI2p — D)~^e, that is, the eigenspaces of A are all of 
dimension 1 and no can be an eigenvalue of A. 

The property (PI) is proved. 

Let us have a look at the localization of the eigenvalues of A. The characteristic poly- 
nomial of A is: 

2p 2p 

xa(a) = n(^^- - ^) - e n (^j- - A) = Q - Q' 

j=l k=lj^k 

w\iemQ = Y\%i{Dj-X). 

In addition to assuming (HI), assume now (H2) is satisfied. Then XA is a real coefficient 
polynomial, so its roots are either real numbers or complex conjugate. Let U be the union 
of the disks of diameters [Dj,D2p-j\, j S i7+- The property (P2) is analogous to Jensen's 
theorem [10]. To prove it, let us consider a complex number z £ C \ U and an integer 
j € J^. Then 

1 ^ 1_ \ ^ _^^^^^mz)-R{D,)r + <^{zr-<^{D,r 



- Dj z-DjJ \z- Dj\'^\z- Dj\^ 

= —Q{z)sj 

where sj > since z ^ U, ^{Dj) being the center of the disk of diameter [Dj,Dj] and Q(Dj) 
its radius. Now assume 2; is a complex, non-real root of XA- Then, taking the imaginary 
part of the equality 

. Q'jz) ^( 1 , i_ 
Q{z) fri\z-Dj z-D,, 

one gets 

= ^(-3(z)s,), 
i=i 

which is absurd. So (P2) stands for all diagonal "conjugate" matrices D. 

Finally, assume (H1)-(H2)-(H3) are satisfied, that is we change D into eD. We also order 
< 0:2 < . . . < ttp-i < OLp. Let A be an eigenvalue of A. If e = 0, A is diagonalizable, its 
eigenvalues are 2p, of multiplicity 1, and 0, of multiplicity 2p—l. According to the theory of 
perturbations, we want to prove that A is necessarily either in a neighborhood of size 0(e) 
of the origin or in a neighborhood of size 0{e) of 2p. Moreover, there should be only one 
eigenvalue, real, in the neighborhood of 2p. We already know that the non-real eigenvalues 
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of A are located in the closed disk P(0, emaxjgj-^(|aj| + Let A be a real root of xa- 

Then a simple computation yields 



1 1 \ „ A — eaj 



= 1. 

One notes at once that necessarily A > eai in order for the sum to be non-negative. Let us 
study the behavior of the function 



y - (^ctj 



Computing the derivative with respect to y, we find that, for e > 0, on the interval 
(e maxjgj-^ (ofj + |/3j|),oo), h{e,-) is decreasing. Since for y > e maxjgj7_^ (aj + |/3j|) > eop, 
h{e, y) > and limoo h{e, •) = 0, there is at most one real eigenvalue larger than eap. Using 
the implicit function theorem, and expanding e i— )• /i(e,A(e)) — 1 in a neighborhood of 0, 
knowing that A(0) = 2p, the only root of XA that is larger than emaxj^j^^aj + |/3j|) can 
be expanded as 

X{e) = 2p (l + ^{a)e + ^{{a- {a)f + P^)e^^ + 0(6^), 
which concludes the proof of (P3). 

□ 

We now turn to the amplification matrix A in (4.1) and express it in terms of P and D. 
This matrix can indeed be written as 



with D = i2peV = e diag(ap — if3p, . . . ,ap + if3p). 

One can then directly formulate the following proposition : 

Proposition 4.2. The eigenvalues A;^^, i = j, . . . , 2p, of the amplification matrix A defined 
in (4.1) are contained in two regions: there are 2p — 1 eigenvalues in the disk 

and one real eigenvalue, an expansion of which is given by 

>^i = l+^-^-^{a' + li'-{af))+St 0(1). 
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The proposition is easily verified by inserting the expression for A into proposition 4.1. 
These eigenvalues can be further examined for the upwind and the centered flux scheme. 
For the centered flux (3.3) for which V = Vc(C)i we have 



so that 



a,- 



0, 



Ax 



6t 



sm\C){v^)+5t o(l), 



5t 
e 



5t 



-Vr 



whereas for the upwind flux (3.2) for which V 



2sin2(C/2) 
2p ' \v 



Ax 



] = ' 

Vc/(C), 

o sin(C) 
= 2p — : V 



Ax 



so that we get, noting that (a) 
5t 



4psin2(C/2)(|t'|)/Ax, 



1 + 



eAx 



2sin^(C/2)(|H) 



Ax2 



sm' 



'(C/2)((t 



sni' 



^(C/2)(|H)2)+<5to 



Ai e ^ I 1 



25t 



j = 2, . . . , 2p. 



We now illustrate this result numerically. We consider equation (2.1) on the velocity 
space (2.2) using p = 10 with e = 1 • 10~^ on a mesh 11 := {xq = — 1 + Ax/2, . . . , 1 — Ax/2} 
with Ax = 0.05 and periodic boundary conditions. We compute the eigenvalues of a forward 
Euler time integration with 5t = and 5t = 0.5e^, respectively, for both the upwind and 
centered flux schemes. The results are shown in figure 3. Clearly, the spectrum of the 
forward Euler time-stepper possesses a spectral gap. The eigenvalues in correspond to 
modes that are quickly damped in the kinetic equation, whereas the eigenvalue close to 1 
corresponds to the slowly decaying modes that survive on long (diffusion) time scales. We 
see that, for both the upwind and central schemes, the fast eigenvalues are centered around 
1 — 5t/e^ . The eigenvalues close to 1 are of order 1 — e for the upwind scheme and of order 



1 



for the central scheme. 



4.3 Projective integration 

The next step is to examine how the parameters of the projective integration method need 
to be chosen to ensure overall stability. It can easily be seen from (3.5) that the projective 
forward Euler method is stable if 



At-{K + l)5t 



At-{K + l)5t 



< 1, 



5t ' y - 

for all eigenvalues \st of the forward Euler time integration of the kinetic equation. 



(4.2) 
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Figure 3: Spectrum of a forward Euler time-stepper for a spatial finite volume formulation of the 
kinetic equation (2.1) for different values of St. Left: upwind scheme. Right: central scheme. The 
inset shows a zoom of the neighbourhood around 1. 

The goal is to take a projective time step At = 0{Ax'^), whereas 5t = O(e^) necessarily 
to ensure stability of the inner brute-force forward Euler integration. Since we are interested 
in the limit e — )■ for fixed Ax, we look at the limiting stability regions as At/6t — )• oo. In 
this regime, it is shown in [8] that the values A^^ for which the condition (4.2) is satisfied 
lie in two separated regions Vf^ U "Cf"^ which each approaches a disk, 

^ V At At J ^ 

The eigenvalues in P^^^ correspond to modes that are quickly damped by the time-stepper, 
whereas the eigenvalues in correspond to slowly decaying modes. The projective inte- 
gration method then allows for accurate integration of the modes in Pf ^ while maintaining 
stability for the modes in ^ . 

Based on the formulae for the eigenvalues A;^^ and the stability regions of projective 
integration, we are able to determine the method parameters St, At and K. The first 
observation is that, to center the fast eigenvalues of the inner time integration (that are in 
P2) around 0, one should choose 6t = e^. Note that this time step is chosen to ensure a 
quick damping of the corresponding modes. The maximal time step that can be taken for 
stability of the inner integration would be 6t « 2e^ due to the bounds in ; in that case, 
however, the fast modes of the kinetic equations are only slowly damped. 

Remark 4.3. For the choice 6t = e^, the spectral properties reveal a natural, but important, 
restriction on the required mesh size Ax, which needs to satisfy Ax > Vp€, to ensure stability 
of the inner forward Euler method. Therefore, the limit Ax — ?• for fixed e is not considered 
in this text. 
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Before deciding on the number of inner forward Euler steps, we need to choose the 
projective, outer step size At. To this end, we require the real eigenvalue of the forward 
Euler time integration to satisfy A^^ G T^i^ , that is 

l_2i^<Ai, <1. 
At ~ ~ 

The second inequality is always satisfied. For the central scheme, we have 

>^st = l-^^ + St o{l). 

Using 5t = and Y^- Vj/{2p) = dp, we then obtain 

At<2-—, 
dp 

which is similar to the CFL condition for a forward Euler time integration of the heat 
equation. (Note that the maximal allowed time step is a factor four larger than that of the 
heat equation.) 

Remark 4.4. The similar derivation for the upwind scheme shows that, in that case. 
At = 0(e), which is undesirable. We will see further on that there are obstructions in the 
consistency analysis too. 

Finally, Ax being fixed beforehand, we need to determine the number of small steps 
K. Introducing r = e/Ax and z/ = dpAt/ Ax'^, stability is ensured if the eigenvalues of the 
forward Euler time integration that are in D are contained in the region Vp. This leads 
to the condition 

VpT < ^ — 

which, after some algebraic manipulation is seen to be equivalent to 

K > 2 \, + ^°^(^/^) 



1 + \og{vp)/ log(r) \og{rVp) ' 

Recalling that Vp = {2p — l)/{2p) and dp = {Ap^ — l)/(12p^), the study of the dependence 
of the bound of -fT in r yields two cases : 

• if < 1/4, that is ra&yip{vp)v < m.mp{dp), then K = 3 independently of r and p, that 
is, if Ax is fixed, independently of e and p; 

• if G [1/4:, 2], if one chooses r < dp/u, K can be safely taken equal to 3 as well. 

Under these hypotheses, we conclude that the projective integration method has an 
e- independent computational cost. 

We illustrate this result numerically. We again consider a forward Euler+centered flux 
formulation of the kinetic equation (2.1) with e = 1 • 10~^ in the velocity space (2.2) using 
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Figure 4: Comparison of the eigenvalues of a forward Euler+centered flux of (2.1) and the stability 
regions of the projective forward Euler method with K = 1 (dashed), K — 2 (dashdotted) and 
A' = 3 (solid). The inset is a zoom on the neighbourhood around 1. 



p = 10 on the spatial mesh IT := {xq = —1 + Ax/2, . . . , 1 — Ax/2} with Ax = 0.05 and 
periodic boundary conditions. The time step is 6t = e^. We plot the eigenvalues together 
with the stability regions of the projective forward Euler method with K = 1,2,3 and 
At = 2Ax^ /dp. The results are shown in figure 4. We see that, in this case, for which 
r = e/ Ax = 0.2, K = inner steps are required to ensure overall stability. Also note that 
the stability region does not depend on K. 



5 Consistency analysis 

Our next goal is to estimate the consistency error in the macroscopic quantity that is 
made at each outer step of the scheme. To this end, we will study the truncation error in Z^, 
keeping in mind the Hilbert expansion f^ = + eg^ (2.3). Again, to simplify notations, the 
subscript e is used only when dealing with the solution of the continuous equation (2.7) 
whereas / is used for the solution of the numerical scheme. Let us recall the brute force 
inner scheme (3.1), in which we now consider vector- valued quantities, hereby omitting the 
subscripts that refer to the dependence in x or v, 

f'^' = Sstf = f' + St (-i^p + , (5.1) 

where $ is the linear spatial discretization operator 

, . _ Hf)i+l/2,j - <P{f)i~l/2,j 

Following the stability condition in Proposition 4.2, we take 6t = and we only consider 
the centered flux (3.3) ; the upwind flux case (3.2) is commented on at the end of the 
section. Note that, with this centered choice, $ satisfies 

('!>(/'=)) = £(<!>(/)). (5.2) 
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The inner scheme then reduces to 

fk+i = S,2f'^ = p^-e^f). (5.3) 
In spatial vector notation = {pf)i^{o,...,Nx}^ the projective scheme (3.5) in p reads 



p^+i = + (At -{K + l)e Y ' ^ , (5.4) 



where + 1 is the number of smah steps, and the values p are obtained as the averages 
over velocity space of the numerical approximations in (5.1). 

Let /e be a smooth solution of (2.7). To bound the truncation error in p of the projective 
integration method (5.4), we introduce the following notations, VA^ > 0, A; € {0, . . . , + 1}, 
(i,j) G{0,...,iVj X J : 

• the exact solution at time t^'^, i.e. /^'^ := fe^j{xi,t^'^)] 

• the corresponding density p^f := p^{xi,t^'^) = (/^^'^); 

• intermediate values obtained through iterations of the inner scheme, starting from the 
exact solution, f^'^ := {8,2^}^; 

• and the corresponding density p^'^ := {f^''')- 

Note that f^'^ = and p^'^ = p^ . The truncation error in p of the projective scheme 
(5.4), is the quantity 

_ pN+l _ pN,K+l ^ / At-{K+l)e^ \ _ pN,K 

At \ At ) £2 • 

To bound in the norm, we first estimate the iterated truncation errors of (5.1) 
with respect to the equation (2.7) for G {0, . . . , K + 1}, given as 

rN,k nN,k 



N,k 
Br :- 



ef'^' = S,2e''^' + -mfj^') - vd^Ut"^^')) + 0{^). (5.5) 



Using a Taylor expansion for the exact solution / combined with equation (2.7), and as- 
suming that / is such that dfif and 9^fe/, k ^ {1, . . . ,K + 1} are bounded uniformly with 
respect to e, a straightforward computation yields 

1 

-I 

e 

We thus have, recalling that e^'*^ = and that the spectrum of S^2 lies in the unit disk, 

^N,K+i ^ 1 ^ gK-k ^^^jN,k^ _ ^5^/^(tA^,^)) + 0{{K + l)e'). (5.6) 

fc=0 

Remark 5.1. In the above formula, as well as in the remainder of the section, to keep the 
notations as clear as possible, the Landau symbol O(-) should be understood as an estimate 
in the L?' norm. 
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Unfortunately, if we compute directly the spatial truncation error in (5.5), the centered 
difference being of order 2, a stiff term Ax^/e appears. We therefore proceed by estimating 
(S'^^^f), where A : /, E C^{T x (0,T);]R2p) ^ _ viOJ^) G M2pxAr, ^^^^^ ~ -g 

again the projection on the discretization points. Note that A is a linear operator and, 
more precisely, that it is the truncation error of the approximation of the first order spatial 
differential by the centered scheme, so that A/^ = 0{Ax'^). In what follows, for the sake of 
simplicity, we will denote the composition # o A (resp. S^2 o A) by the product ^>A (resp. 
S,2A). 

The crucial first step is to see that (5.3) reads 

5,2A/, = (A/,) - ea>A/„ 

and, consequently, 

A/, = (A/,) - e{(<^A/,) + 1>(A/,)} + e^<^^AU 

A simple combinatoric argument implies that, for A; > 3, 

S'^^Af, = (A/,) - e{{<^Af) + $(A/,)} 

+ {(cD^A/,) + <l>^{Af) + <I>(#A/) + {k- 3)(«>2(A/))} + 0{e^). 

Taking the mean value of S^2Af^, and using (5.2), as well as the fact that the linear operator 
$ — vdx is odd in v, we get : 

• for A: = : (A/,) = e{Ag,) = eO{Ax^), 

• for A: = 1 : (5,2 A/,) =e{{Ag,) - ($Ap,)} - e^{<^Ag,) = eO{Ax^) + e^O{Ax^), 

• for /c > 2 : 

{S^,Af,) =e{{Ag,) - (4>Ap,)} + {{<l>^Ap,) - {^Ag,)} + 0{e'), 
= eO{Ax^) + e^O{Ax^) + 0{e^). 

Using Young's inequality and plugging this estimate in (5.6), we get 
and, consequently. 



At \ At J € 

At-{K + l)e^ \ (e^'^+i - e^'^) 



c2 



At 

At J ' y ' ^' At) 

At 



\-{K + X)^\ d,p{t^'^^') -(1- {K + d,p{t^^^^^) 



+ 0{At) + e^O ('^4^ ) + 0(6^) + 0{Ax^). 
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Following the classical parabolic CFL condition At = O(Ax^), we get 



E^ = 0{At) + 0(^]+0{e'). 




The projective scheme is consistent with (2.5) at order 2 in e and, as e goes to 0, the limiting 
scheme is consistent at order 1 in time and 2 in space with (1.2). 

Remark 5.2 (Upwind fluxes). For the upwind flux, the fact that ($(pe)) = O(Ax^), instead 
of vanishing as in the centered flux case, implies a consistency error term of order Ax^/e 
that is not easily cancelled. 

Remark 5.3 (Hilbert expansion). When taking 6t = e^, rewriting the scheme (5.1) in 
terms of p and g = {f — p)/e leads to 



which is a scheme for (2.4). From this equation, the effect of the particular choice of time 
step 5t = becomes clear. With this time step, in accordance with (2.3) and (2.6), we see 
that the distribution satisfies 



and, therefore, the projective integration scheme recovers the first two terms in the Hilbert 
expansion of /. 

In conclusion, we summarize the above results on stability, consistency and the number 
of steps : 

Theorem 5.4. Under the CFL condition At = Ax^ /{Adp), for K > 3, assuming d^^f and 
Q^kf > k S {1, . . . , A' + 1} , are bounded uniformly with respect to e, the following estimate 
holds for all NAt < T, 



The estimate (5.8) shows that, if e is fixed, the optimal choice of At in terms of accuracy 
is At = O(e^). Of course, this leads to prohibitive costs as e — )■ 0, but it allows us to consider 
a solute computed with this choice of time-step to be precise at order e^. Larger values of 
At, and in particular the values At = 0{Ax'^) that we envision, increase the error due 
to the outer step; smaller values increase the error due to the time derivative estimation 
from the inner steps. Note also that taking At — t- for fixed e would completely defeat 
the purpose of the projective integration method. Additionally, in this limit the method is 
unstable, since this choice implies Ax — t- for fixed e (see also remark 4.3). 



Pi 



k+l 



p\ - e2(c^,(/)) 

-$.(p'=)- 6 {$,,(/) -($,(/))} 




= + eg"" =p^ - €^{p^) + 0(e2) 




(5.8) 
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6 Numerical results 



In this section, we illustrate the above results via the numerical simulation of two model 
problems, namely the linear kinetic equation (2.1), and the Su-Olson problem (2.8). 

6.1 Linear kinetic equation 

We consider equation (2.1) on the velocity space (2.2) using p = 10 with e = 5 • 10~^ on the 
spatial mesh IT := {xq = —1 + Ax/2, . . . , 1 — Ax/2} with Ax = 0.1 and periodic boundary 
conditions. As an initial condition, we take 

{2, for - 0.5 < X < 0.5 and - 0.75 <v< 0.25, 
_ _ _ _ 

1, otherwise. 

We define the initial condition for the recursion (3.1) by taking cell averages, i.e. 

fi,j= / Mx,v,t = 0)dx. 

Jxi-Ax/2 

We perform a time integration up to time t = 2.5 using a centered flux/forward Euler scheme 
with 6t = e^, as well as a projective forward Euler integration, again using 6t = e^, and 
additionally specifying K = 4 and At = vAx /dp, with dp = {v^) = 0.3325 (using the cell 
averages) and i/ = 1. For comparison purposes, we also compute the result using a centered 
flux/forward Euler scheme with 5t = e^, which we will consider to be the "exact" solution, 
and a solution of the limiting heat equation on the same mesh using At = OAAx'^ /dp. The 
experiment is repeated for e = 2 • 10~^ and Ax = 0.1 and Ax = 0.05, respectively. The 
results are shown in figure 5. We show both Pe{x,t = 2.5) and the flux J^{x,t = 2.5), with 

Je(x,t) = - / vfe{x,V,t)dv. 

e Jv 

We see that the complete simulation with 6t = and dt = visually coincide in all cases. 
The projective integration method results in a solution that is closer to that of the limiting 
heat equation. Note that the differences between all solutions become smaller for decreasing 
e and Ax. (For Ax = 0.1, the difference between projective integration and the "exact" 
solution are mainly due to the space, and correspondingly large time step, in accordance 
with theorem 5.4.) 

Next, we look at the convergence properties in terms of e. To this end, we repeat the 
computation for Ax = 0.1 and different values of e. We consider the centered flux/forward 
Euler flux finite volume scheme with dt = to be the "exact" solution and approximate 
the error of the other simulations by the difference with respect to this reference simulation. 
(This reference solution is at least an order in e more accurate than the other results.) We 
first investigate the error of the centered flux/forward Euler flux finite volume scheme with 
St = e^. The results are shown in figure 6 (top). The O(e^) behaviour is apparent. In the 
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Figure 5: Results of simulation of the kinetic equation (2.1) at time t = 2.5. Left: density; right: 
flux. Parameter values are (a) Aa; = 0.1 and e = 0.05 (top); (b) Aa; = 0.1 and e = 0.02 (middle); and 
(c) Aa; = 0.05 and e = 0.02. Shown are (i) a centered flux/forward Euler flux scheme with St = 
(solid); (ii) a centered flux/forward Euler scheme with St ~ (dashed); and (iii) the projective 
integration method (dashdot). For comparison, also the solution of the heat equation is shown 
(dot). 
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Figure 6; Top: Error in density (left) and flux (right) of a centered flux/forward Euler flux integration 
of the kinetic equation (2.1) using Aa; = 0.1 and St = as a function of e at time t = 1.25 (dot), 
t = 2.5 (dash) and t = 3.75 (sohd). For comparison, we also plot a line with slope 2 (dashdot), 
indicating the error that is predicted by the consistency analysis. Bottom: Error in density (left) and 
flux (right) of a projective forward Euler integration of the kinetic equation (2.1) using Aa; = 0.1, 
St ~ e^, K = 3 and At = Ax'^/dp as a function of e at time t — 1.25 (dot), t = 2.5 (dash) and 
t — 3.75 (solid). The error is 0{Ax^) was predicted by theorem 5.4. 
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Figure 7: Error in density (left) and flux (right) of a projective forward Euler integration of the 
kinetic equation (2.1) with e = 2 • 10"^ using Ax — 0.1, St — e^, A' = 3 as a function of At at time 
t = 1.25 (dot), t = 2.5 (dash) and t = 3.75 (solid). For comparison, we also plot a line with slope 1 
(dashdot), indicating the error that is predicted by the consistency analysis. 

same way, we consider the error of projective integration using K = 3 and At = Ax'^/dp. 
Figure 6 shows that this error is largely independent of e, especially when e — t- 0. 

Finally, we look at the error of projective forward Euler as a function of At. We perform 
a projective forward Euler simulation using Ax = 0.1 and e = 2 • 10~^ using 5t = and 
K = 3. As the projective step size, we use At = uAx^ /dp for a range of values of i^, and we 
again compare the density and flux with respect to the reference solution. The results are 
shown in figure 7. We clearly see the first order behaviour as At — >• 0. Also remark that, 
for large values of At, the error increases quickly due to a loss of stability. From this figure, 
it can be checked that the loss of stability indeed occurs at = 2, as discussed in section 
4.3, whereas for the limiting heat equation, v = 0.5 is the maximal value that should be 
used to ensure stability. 

6.2 Su-Olson problem 

We consider equation (2.8) on the velocity space (2.2) using p = 10 with e = 5 • 10~^ on 
the spatial mesh 11 := {xq = — 1 + Ax/2, . . . , 30 — Ax/2} with Ax = 0.1 and homogeneous 
Neumann boundary conditions. As an initial condition, we take 

f,{x,v,t = Q) = e,{x,t = Q)=A, 

with A = 1 and A = 1 • lO"^'^, respectively. Again, we take the cell averages of the initial 
condition, and perform a time integration up to time t = 1 using (i) a centered flux/forward 
Euler flux scheme with 6t = e^, and (ii) a projective forward Euler integration, again using 
5t = e^, and additionally specifying K = 3 and At = uAx'^/dp, with dp = {v^) = 0.3325 



22 




0.01 0.1 1 10 100 0.01 0.1 1 10 100 0.01 0.1 1 10 100 



Figure 8: Results of simulation of the kinetic equation (2.8) at time t = 1. using A = 1- 10^^" (top) 
and A = 1 (bottom). Left: Pe{x,t ~ 1); middle: 6e{x,t ~ 1); right: eJe{x,t — l)/p^{x,t ~ 1). Shown 
are (i) a centered flux/forward Euler flux scheme with St — (solid); (ii) a centered flux/forward 
Euler scheme with St = (dashed); and (iii) the projective integration method (dashdot). All 
simulations visually coincide. 

(using the cell-averaging) and v = 1. For comparison purposes, we also compute a reference 
solution using a centered flux/forward Euler flux scheme with 5t = e"^. The results are 
shown in figure 8. We show Pe{x,t = 1) and 9e{x,t = 1), as well as the quantity eJe/pe-, 
which is a measure of the "limited-flux property", that is, if is nonnegative, we should 
have 

e\Je 

Only the result of projective integration is visible, since the solutions using the other pro- 
cedures are visually indistinguishable on this scale. The errors of the full forward Euler 
simulation and the projective integration are shown in figure 9. We clearly see that, while 
the computational cost of projective integration is much lower, the error remains of the 
same order of magnitude. 

7 Conclusions and discussion 

We investigated a projective integration scheme for the numerical solution of a kinetic equa- 
tion in the limit of small mean free path, in which the kinetic description approaches a diffu- 
sion equation. The scheme first takes a few small inner steps with a simple, explicit method, 
such as a centered flux/forward Euler finite volume scheme, and subsequently extrapolates 
the results forward in time over a large outer time step on the diffusion time scale. We pro- 
vided a stability and consistency result, showing that the method is asymptotic-preserving. 
We conclude with some remarks, and some directions for future results. First, a higher 
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Figure 9: Error of the numerical solution of the kinetic equation (2.8) at time t = 1. using A = 
1 ■ 10-1° (top) and A = 1 (bottom). Left: p^{x,t = 1); middle: e^{x,t = 1); right: ej^{x,t = 
1)/ Pf:{x, t = 1). Shown are the error with respect to a reference solution of (i) a centered flux/forward 
Euler flux scheme with 5t ~ (solid); and (ii) the projective integration method (dashdot). 



order outer integration method can readily be used to obtain a higher order in the macro- 
scopic time step At, see e.g. [23, 19]. We emphasize that this higher order accuracy does 
not depend on the order of the inner simulation, since the error in time at that level is of 
the order 0{e^), due to the choice of the inner time step. 

Several new directions are currently being pursued. First, we are extending these results 
to the kinetic Fokker-Planck case, which requires a precise study of the discretization of the 
second-order derivation in the velocity variable. Furthermore, we are also considering pro- 
jective integration in conjunction with a relaxation method [1] to obtain a general method 
for nonlinear systems of hyperbolic conservation laws in multiple dimensions. 
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